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Given for one instant an intelligence which could comporehend all the forces 
by which nature is animated and the respective situation of the beings who com- 
pose it-an intelligence sufficiently vast to submit these data to analysis-it would 
embrace in the same formula the movements of the greatest bodies of the universe 
and those of the lightest atom; for it, nothing would be uncertain and the future, 
as the past, would be present to its eyes.-Laplace. 



I. INTRODUCTION 



Simulating chemical reactions in the condensed phase involves the parallel computation 
of the concerted quantum dynamics of breaking and formation of chemical bonds with 
the simultaneous response and influence of the surrounding media. In many systems it is 
not entirely possible to treat all of the motions quantum mechanically without resorting 
to various levels of approximation. One approach which has proven to be quite powerful 
in a wide range of scenarios is to partition the full system into interacting subsystems 
according to whether the dynamics in the subsystem is predominantly classical-like or if 
the dynamics must necessarily be treated quantum mechanically. Such partitionings are 
motivated by the simple observation that quantum mechanical effects are typically localized 
to only a few specific types of motions typically along some chosen reaction coordinate. 
The remaining motions, i.e. those of the surrounding media, can largely be regarded as 
classical. Thus, partitioning the extended system into interacting quantum and classical 
subsystems allows one to study the details of a condensed phase quantum phenomena while 
retaining a detailed description of the response and influence of a complicated molecular 
media in its full dimensionality. Consequently, the level of description afforded by mixed 



quantum/classical treatments and the relative computational efficiency of the method has 
allowed a synergistic interaction between theory and experiment since experiments can be use 
to validate theoretical models which in turn have provided a number of valuable mechanistic 
insights. 

Current quantum chemical methods have proven to be quite adequate in providing struc- 
tural and dynamical information for radicals in the gas phase and inert matrices. In these 
circumstances, quantum computations have provided a direct link between experimental 
spectroscopic investigations limited chiefly by the level of sophistication needed to obtain 
accurate structures, spin dependent properties such as hyperfine interactions, and higher 
order propertied such as hyperpolarizabilities. 

In the condensed phase, the story is quite different. In order to perform molecular dy- 
namics of the solvent molecules, we need to provide accurate empirical solvent-solvent and 
solute-solute interaction potentials However, empirical models cannot accurately re- 

produce many important modifications induced by the solute on the electronic structure 
of the solvent molecules. Consequently, the use of simple fixed point charge models for the 
solvent ignores the polarization response of the solvent molecules to changes in the electronic 
structure of the solute. Mixed approaches have been developed, in which a Monte Carlo 
simulation of the solvent was coupled to a quantum mechanical description of the solute, 
however, to date these have been performed at the Hartree-Fock level due the computational 
overhead required to sample a statistically meaningful number of solvent configurations . 
Furthermore, the representation of the solvent molecules by simple point charge models ne- 
glects many important polarization induced effects when computing the electronic structure 
of the solute. 

Solvent models have been developed in which the solvent is treated as a polarizable di- 
electric continuum 1^-0 . While the use of dielectric continuum models has a long history, the 
latest generation of models are extremely sophisticated and include such features as realistic 
solute cavities and exact description of the electrostatic problem Such Polarizable 

Continuum Models (PCM) have been included in the Gaussian/94 (G94) |jT2| package al- 
lowing an accurate extension of a variety of quantum chemical approaches (including HF, 
MP2, CC, and DFT) to condensed phase problems. Such fully ab initio approaches are 
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currently accurate enough to interpret experimental hyperfine splittings and permit an un- 
biased assessment of the role of electronic, vibrational, and solvation effects in determining 
the observed result |10|- The success of these models open the possibility for studying a 
number of large biologically important radical species in their natural aqueous environment. 

Ab initio molecular dynamics methods, especially the Density Functional Theory- 
Molecular Dynamics (DFT-MD) methods introduced by Carr and Parrinello [|l3|-[l^ (CP) 
have been highly successful in describing a variety of chemical phenomena in the condensed 
phase [P!3|-T7|. The CP method avoids the necessity of providing semi-empirical pseudopoten- 
tials for the molecular forces by solving the Kohn-Sham equations for the electronic density 
for a given nuclear configuration. The molecular dynamics are then solved under the Born- 
Oppenheimer approximation in which the nuclei move under the Hellmann-Feynman forces 
computed by the DFT calculation. Additional quantum effects in the nuclear motion, such 
as energy level quantization and barrier tunneling, can be introduced through path integral 
dynamics |TB|,|T^. While this methods is useful in determining mechanistic details of many 
chemical reactions, electronic transition due to the non-adiabatic motion of the nuclei can- 
not be properly include due to the fact the Kohn-Sham equations are not guaranteed to be 
convergent variationally for states other than the electronic ground state. Thus, all nuclear 
dynamics are confined to the ground adiabatic potential energy surface. For a large number 
of systems of chemical importance, this level of description is sufficient. However, for many 
systems, sudden changes in the solute quantum state play pivotal roles in determining the 
dynamics of the reaction. 

Quantum-classical simulations overcome many fundamental shortcomings associated 
with both classical molecular dynamics (MD) and the quantum chemical methods men- 
tioned above. Foremost amongst the assumptions in both MD and quantum chemistry is 
that all the dynamics occur on a single potential energy surface under the Born-Oppenheimer 
approximation. Quantum mechanical effects arising from tunneling and dynamical interfer- 
ence effects and discrete energy levels are completely absent in classical MD. These effects 
become significant particularly in low temperature systems containing hydrogen atom mo- 
tions and for computing various time-dependent spectroscopic observables. The restriction 
of the dynamics to a single Born-Oppenheimer potential surface is another serious limita- 
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tion. For each non-degenerate electronic state of the many-atom system there is a distinct 
Born-Oppenheimer potential energy surface parameterized by the positions of all the atoms. 
When a transition between electronic states occurs, the forces acting on the atoms may 
change dramatically. Thus the proper and self consistent incorporation of these effects is 
crucial to describing a wide host of chemical dynamical effects, many of which are important 
to free radical dynamics. ||20[| 

Electronically non-adiabatic processes play a surprisingly ubiquitous role in chemical 
reaction dynamics for a wide number of systems including photosynthesis, polymeric chain- 
reactions, photodissociation, and electron transfer reactions, and various aspects which 
must be considered in order to simulate such processes for a condensed matter system. 
Foremost amongst the goals of this chapter is to present the quantum-classical methodology 
as a useful method of simulating these processes in the condensed phase, to discuss some of 
the subtle aspects of non-adiabatic dynamics which we have demonstrated to be important 
considerations, and to present some of our recent results obtained using these methods. 

We consider here a system composed of a few important degrees of freedom which include 
any reaction coordinate we care to choose coupled to a much larger number of degrees of 
freedom which we shall refer to as the bath. We shall treat the motions of the system using 
quantum mechanics and motions in the bath via classical mechanics. The quantum motions 
may be, for example, nuclear motions, as in the hydrogen transfer reaction, electronic de- 
grees of freedom as in the case of an electron transfer reaction or non-radiative relaxation 
process, a combination of the two, and so forth. One important class of motions of chemical 
import which are inherently quantum mechanical in nature are transitions between discrete 
electronic states. Since the nuclear positions and velocities of all the molecules in the system 
determine the electronic states and the transition amplitudes between the states and the 
states and transition amplitudes in turn influence the the forces governing the trajectories, 
any mixed quantum classical theory must treat self-consistently the classical and quantal 
degrees of freedom. 

The aqueous electron, nature's simplest radical anion, provides an ideal test case system 
for the methods which we shall develop below. Many of the dynamical features present in 
this system have analogues in other systems in which energy quantization, non-radiative 
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energy transfer, and other quantum effects play import roles. The aqueous electron is one 
of the very few important chemical systems which can be described more or less exactly 
without resorting to various levels of quantum chemical approximations. Thus, it provides 
a unique testcase for a variety of dynamical theories. 

The remainder of this chapter is organized as follows: we first present a general overview 
of classical molecular dynamics and Monte-Carlo methods. This is aimed as a general 
overview of the use of molecular simulational methods. As a part of this overview, we 
discuss recent developments in ab initio molecular dynamics, in particular the density func- 
tional theory methods developed by Carr and Parrinello. We then focus our discussion on 
nonadiabatic simulations. We present the formal details of quantum/classical dynamics and 
highlight two important considerations which are necessary to address: the treatment of in- 
teraction force and the treatment of the phase coherence between the quantum and classical 
subsystems. We then move on to show results of simulations based upon these methods fo- 
cusing primarily upon the role that quantum decoherence plays in determining the survival 
timescale of an excited electronic state using the aqueous electron as a prototypical system, 
highlighting the interplay which has evolved between theoretical and experimental studies 
of this system over the past few years. Finally, we conclude with a discussion of future 
directions we and others are undertaking to apply the methods discussed here to a variety 
of other systems. 



II. EXCESS ELECTRON AS A PROTOTYPICAL FREE RADICAL SYSTEM 



The solvated electron provides perhaps the simplest example of electronic dynamics in 
the condensed phase. Aqueous electrons are of particular importance since they play a 
prominent role in the broad area of radiation chemistry of water. In spite of intensive 
experimental and theoretical study since identification of the hydrated electron over 30 years 
ago, many dynamical features of electron solvation have remained incompletely understood, 
largely due to the extremely fast timescales on which electronic relaxation occurs ^l|] . Newly 
developed femtosecond spectroscopic techniques are now providing glimpses of the details of 
condensed phase electronic dynamics. From a theoretical standpoint, the aqueous electron is 
an ideal test case for new condensed phase theories. Recent theoretical advances in treating 
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condensed phase electronic dynamics p^-p6| has led to an exceptional interplay between 
theory and experiment, helping to unravel many subtle dynamical aspects of this system. 

While excess aqueous electrons are often prepared in conjunction with chemical reactions 
and charge transfer reactions in H2O, experimental preparation of the electron in neat water 
can be accomplished by direct photoionization of water. Following injection into the water, 
the electron typically has 2.0 to 2.5 eV of excess energy which must be absorbed by the 
water as the electron is solvated. Initially, the electron is delocalized in a continuum band of 
energy levels. However, the electron rapidly localizes into one of three nearly degenerate p- 
state and subsequently relaxes to an s-state where it is further solvated by the surrounding 
water molecules. This fully solvated species exists until scavenged by impurities in the 
system (typically on the order of a few hundred ns for ultra-pure water). The kinetics of 
the localization process has been extensively studied both theoretically and experimentally. 

The fully solvated electron can be regarded as a particle in a hard roughly prolate 
spheroidal cavity formed by the surrounding water molecules. The first "solvation shell" 
of water molecules is an octahedral structure with hydrogen atoms on the vertices with 
an average e-/0 radius of ~ 3.2 A. In the fully solvated state, the water molecules in the 
surrounding shells arrange themselves as to maximize the hydrogen bonding network. The 
first excited states of this system are p-states lying roughly 1.7 eV above the ground state 
which are split due to slight asymmetry of the shell of water molecules about the electron. 
The energy gap between the s and p states and the energies of the p states are modulated 
by fiuctuations in the local solvent structure. 

In the simulations which we shall discuss through the course of this chapter, we assume 
that the electron has been prepared in a equilibrated s-state and subsequently promoted 
into one of the p-states when the s-p energy gap becomes 1.7 eV simulating the effect of 
the 720 nm laser excitation pulse used in the experiments performed by the Barbara and 
co-workers. The relaxation mechanism following excitation is given by 

hu ^ tnr * /-| \ 

Following promotion into one of the p-states at t = 0, the system is no longer in equilibrium 
and solvent water molecules move to accommodate the new distribution of electronic charge 
and partially solvates the p-state on a timescale, ti. During this time, the s-p energy gap to 



narrows from 1.7 eV to uJsp ~ 0.5 eV. As the energy gap narrows, the non-radiative transition 
probabihty back to the ground state increases and is maximized around avoided crossings. 
Thus, measurements of the excited state hfetime reflect both the non-radiative decay of the 
state as well as the solvation/relaxation dynamics of the surrounding media. The solvation 
dynamics of this system have been extensively studied by Schwartz and Rossky. The excited 
state solvation response, defined as the normalized autocorrelation function of the s-p energy 
gap following excitation, 

^(t) = / "-f\-"-^^°\ \, (2) 
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indicates a rapid ~ 24 fs inertial Gaussian component which makes up roughly 40% of the 
response followed by a much slower 240 fs exponential decay which makes up the remainder. 
Following the non-radiative p — s transition, the new s-state (e*) is no longer in equilibrium 
and the water solvent molecules respond to resolvate the electron. Figure |1] shows the results 
of two typical runs following preparation in the excited state. 

An important puzzle concerning the non-adiabatic dynamics of the hydrated electron 
has been the surprising lack of a sizable isotope effect on the non-adiabatic transition rate. 



| 27| , |28[| Estimating the non-radiative coupling between the excited and ground state using 



Fermi's Golden Rule 

kio oc \v.{iJo\'^i'ex)\'^ (3) 

shows that the nuclear velocities, v, play a direct role in determination of the non-adiabatic 
transition rate. Since the fastest nuclear velocities in D2O (the D's) are classically \/2 times 
slower than those in H2O (the H's) while the other factors (the electron-water interaction 
potential, the quantum force on the nuclei, etc.) remain the same, the expectation is that 
non-adiabatic transition rates should be roughly half as large in D2O compared to H2O. 
Indeed, mixed quantum-classical simulations have suggested isotope effects of factors of 2 — 



4 for the electronic transition rate in this system. [29-31[ Experiments, however, have found 



at most a modest difference in the non-adiabatic transition rate for electrons photo-injected 



into H2O versus D2O, |32-34| and at most a 15% isotopic effect have been extracted from 
femtosecond spectroscopic studies of photoexcited equilibrium electrons. p8| , p5| This dispar- 
ity between theory and experiment comes as a surprise given the level of sophistication of the 



theories and their abihty to match experiment on many other aspects of the system. This 
has caused us to re-evaluate the quantum-classical methodology and examine the various 
approximations used to derive the method, which we shall describe in detail in the following 
section. 



III. QUANTUM-CLASSICAL METHODOLOGY 

A. Non-adiabatic Dynamics 

Due to the fundamental role that dynamical process play in chemical reactivity, there 
is a large body of literature devoted to extracting dynamical information from computer 
simulations for both adiabatic and non-adiabatic dynamics. [^,^,^^44| As discussed 



above, it is not computationally feasible to treat a condensed phase dynamical system fully 
quantum mechanically while retaining a molecular level description of the dynamics; it 
is, however, entirely possible to simulate the physics of interest by treating a select few 
degrees of freedom quantum mechanically while treating the remainder classically. The 
general prescription goes as follows: write the total Hamiltonian in terms of "quantum" and 
"classical" contributions,^ 

Ht,, = HQ{x,R) + ^ + Vc{R), (4) 

where Hq{x,R) is the Hamiltonian for just the quantum variables, x, parameterized by 
the classical coordinates, R. The other are the classical kinetic energy and the solvent- 
solvent potential interaction, Vc{R), which is independent of the quantum subsystem. The 
quantum subsystem is propagated in concert with the classical particles using the time 
dependent Schrodinger equation, (note: we set h = 1 throughout) 



^Rigorously, one should start from fully quantum realization of the system and introduce mixing 
between the electronic and vibrational motions through the application of the vibrational kinetic 
energy operator on the electronic wavefunction to arrive at these equations. C. f. M. Born and K. 
Haung, Dynamical Theory of Crystal Lattices, (Oxford Univ. Press, 1954). 
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dt ^ "^dR^ 



^{x;R,t) = HQ{x,R)^lj{x;R,t), (5) 



where the summation is over all classical coordinates. Here we have used the chain rule 
of differentiation to decompose the the total time differential, to show explicitly the "non- 
adiabatic coupling term" : RV Rip{x\ R, t). Under the Born-Oppenheimer approximation, the 
quantum subsystem is assumed to respond much faster than the timescale set by the motion 
of the classical particles. When this is the case, the non-adiabatic coupling term can be 
ignored and the classical particles evolve on the adiabatic potential energy surfaces defined 
by solving the time independent Schrodinger equation at each classical configuration, R, for 
the adiabatic state 

HQiR,t)\MRm=eniR,t)\MRm- (6) 
Under the adiabatic approximation, the classical particles evolve under 

mMt) = -^(^"(^W) + Mm))- (7) 

Thus, the coupling between solvent and solute are those predicted by the Hellmann-Feynman 
theorem 

F^m)) = -^{MRmHQiRmMRim 

When the non-adiabatic coupling becomes significant, such as at an avoided crossing of 
the adiabatic energy surfaces, the classical motion can induce transitions in the quantum sub- 
system. As the classical particles move away from regions of strong non-adiabatic coupling, 
the classical particles must evolve on a given adiabatic surface. This asymptotic condition 
introduces a twist on the classical dynamics since the classical particles must "switch" sur- 
faces during the course of a quantum transition. Various computational schemes have been 
developed which incorporate this switching aspect into the classical dynamics. Perhaps the 



most straightforward scheme is the "surface hopping algorithm" pioneered by TuUy |38,40 
subsequently incorporated in the work of Webster, Friesner, and Rossky [^3|,^ and of Coker 
and coworkers |]37| , ^ , |45| , ^ . The general feature of all these algorithms is that the classical 
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particles stochastically "switch" between quantum states according to a criteria determined 
by the quantum transition probability between the initial and final quantum states over a 
given time interval. Thus, over the course of a simulation, a switching trajectory emerges as 
the quantum system hops between the various quantum states tracing out a path through 
state space. 

In mixed quantum-classical simulations, we assume that the nuclear motion of at least 
a few degrees of freedom can be treated via classical mechanics and that the evolution of 
remaining degrees of freedom can be treated via quantum mechanics parameterized by the 
classical variables. The coherent evolution of the reduced density matrix along a given 
classical trajectory, p{R), is given by the Liouville equation 

i^^^^[HQiRit)),piRm- (9) 
Here, the density matrix, p{R{t)), is defined along a single trajectory, R{t) as 

p{R{t)) = \iJ{Rm{ij{Rm, (10) 

where \tp) is the solution of the time dependent Schrodinger equation. In the basis of 

adiabatic eigenstates this becomes 



p{R{t)) = j2p,,{RmMRm{um)i (n) 

with the time evolution of each component given by 

iPij — {^i^jk — ^Rdik) Pkj — Pik {^Ak — ^RiJ.dkj) ■ (12) 

where are the adiabatic energies computed for a given point along the path, dfj is a 
component of the non-adiabatic coupling vector 

d^j^iMRim-^jmm, (13) 

and {\(l)j{R{t)))} are the adiabatic eigenstates of the quantum Hamiltonian at configuration 
R(t). 

Over a short time interval 6t, we can construct the probability for making a quantum 
mechanical transition from an initial state \4>i{t)) to a final state \4>f{t + St)) writing 
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rto+St 

Pij{to + St, to) = \ {(f)j \ exp(-i / dsH{s))(p^ 



Pii[t) 

« <5t^4^ (14) 

The first line is the definition of a quantum transition probabihty The second and third lines 
are obtained by solving the equations of motion for the density matrix over a short time 
period which gives the probability of switching from state i to j by the rate of change of the 
population of the final state weighted by the population of the initial state integrated over 
the short time interval. Using this we can write 

Pii\^) 

These probabilities are computed at each time step over the course of the simulation and de- 
termine whether or not the classical degrees of freedom will continue to evolve on the present 
energy surface or if they must switch energy surfaces due to a non-adiabatic transition in the 
quantum degrees of freedom. This particular approach for computing the quantum transi- 
tion probability has the desired advantage that it minimizes the number of switches between 
states all the while maintaining the correct statistical distribution of state populations at 



all times. |^ 



In the Molecular Dynamics with Quantum Transitions (MDQT) algorithm P(]|J^ , the 
classical variables switch suddenly between the adiabatic potential surfaces at the switching 
times. Such abrupt changes in the potential energy when an electronic transition does occur 
causes a change in the total energy of the system. In the MDQT method, the velocity of 
the classical variables are modified in order to properly conserve energy by randomly 
rescaling the velocity vectors of the classical particles along the direction of non-adiabatic 
transition vector such that the transition energy is transferred from the quantum to the 
classical degrees of freedom (or vise verse) in such a way that the total energy of the system 
is conserved. We shall next discuss a theory in which the trajectories move on an effective 
potential which changes smoothly as the quantum state evolves from the initial to final state 
over a finite time interval. 
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1. Solvent- Solute Coupling in the Presence of Quantum Transitions 



When a quantum transition does occur, the forces couphng the quantum system to the 
solvent change to reflect the new electronic state. This change is often dramatic, especially 
in transitions between an excited state and the ground electronic state. For example, in 
the case of an excess electron in H2O, the first excited state is a p-state while the ground 
state is a roughly spherical s-state. Furthermore the partially solvated p-state occupies a 
much larger solvation cavity than the ground state for the same solvent configuration. Thus, 
both the shape and spatial distribution of the electron changes suddenly as the result of an 
electronic transition. Consequently, the forces coupling the electron to the solvent change 
suddenly as well. A proper description of the solvent-solute coupling must accommodate 
these sudden changes. 

A further consideration is the fact that the total energy of the system must be conserved. 
Any energy involved in making a transition in the quantum state must be absorbed or 
transferred to the solvent continuously throughout the transition. In the case of an excess 
e~in H2O, the p-s transition dumps roughly 0.5 cV of electronic energy into the solvent. 
This energy must be accommodated by the solvent translational, librational, and vibrational 
motions in a time scale of a few femtoseconds. 

Along a given path, R{t), the partial quantum transition amplitude between an initial 
and final quantum state, \i{Ri)) and \ f{Rf)) respectively, is given by 

T,f[R(t)] = (/(it:/)|e-'^''^'''^'''^^^l|i(it:,))e+'^[^(*)l 

= %[i?(t)]e+^^[^(*)l, (16) 

where Uif[R{t)] is the transition amplitude for the quantum subsystem and Sc[R{t)] is the 
classical action computed along the path R{t), i.e. 

Sc[R{t)] = dt[^R{tr - Vc{R{t))] . (17) 

We are calling this a partial amplitude here in the sense that this represents the full transition 
amplitude along one Feynman path taken by the solvent or bath degrees of freedom. Indeed, 
the full quantum transition probability of starting in some initial solvent/solute state \i{to)) 
ending up in a given final solvent/solute quantum state, \f{tf)), at some time tf is computed 
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by summing over all possible Feynman paths, R{t), taken by the solvent between initial and 
final solvent configurations, Ri and Rf, and with initial and final electronic states, \i{Ri)) 
and \f{Rf)) respectively. I.e. 

Kj,{Rj,R,-t) = J VR{t)T,f[R{t)] 

= I Di?(t)(/(/?;)|e-^^''^'''^f''^^^'|i(i?.))e+^^^[«W] (18) 

The transition amplitude for the quantum subsystem, Uif[R{t)], is a complex functional 
of R{t). We can obtain a semiclassical approximation for the combined system/bath propa- 
gator, Kfi{Rf, Ri] t) on the assumption that the magnitude of U[R\ varies much more slowly 
than its phase along a given path. Expanding the phase to to second order in the classical 
paths R yields the variational equation 

5 (hIm\nUif{R) + Sc{R)) = 0. (19) 

The variation of the electronic propagator is given by 

5Uif[R] = -y^'' ds{f{tf)\U{tj,s)SHQU{s,UMU)). (20) 

Using this, we arrive at the "classical" equations of motion for the saddle point trajectories, 
R. 

(We shall drop the tilde for now on.) This is basically Newton's equations of motion. Thus, 
the semiclassical estimate for the force contribution between the quantum and classical 
degrees of freedom is given by 



This force is a functional of the trajectory which we are trying to calculate over the course 
of a simulation and must be solved iteratively. In other words, we guess a path and a set 
of end points, compute the quantum wavefunction for the new configuration, propagate 
the initial wave forward in time and the final wave back in time to a set of intermediate 
times to compute the forces for a new path and continue with this cycle until either the 
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procedure satisfies some convergence criteria or one simply runs out of patience or CPU 
time. Tliis semiclassical/stationary pliase description of the coupling between the quantum 
and classical variables is a central feature in the quantum/classical stationary phase surface 
hopping methods developed by Webster, Friesner, and Rossky (WFR) [^,^,^]. While 
providing an exact semi-classical description of the forces coupling the quantum and classical 
subsystems, the computational overhead required to iteratively determine the classical path 
can be quite large, especially if there are two or more paths with nearly identical action or if 
the time interval is taken to be too large. This creates considerable computational overhead 
since the motion of all of the solvent species must be determined variationally rather than 
through a propagative scheme, such as the Verlet algorithm. 

Furthermore, when we specify the initial and final states of the quantum solute, we are 
in effect imposing a quantum measurement on the solute which destroys coherence in the 
quantum subsystem. In the algorithm developed by Webster, Friesner, and Rossky (WFR), 
this "computational coherence timescale" is chosen to be identical to the timescale required 
to accurately compute the fastest motions in the solvent. Later developments by Coker and 



co-workers ||37| , ^ , ^ , ^ and later by Bittner and Rossky ||36| , ^8| -pO[] have allowed for full 
retention and partial retention of coherence. 



2. Treating the phase coherence between the solvent and solute subsystems: Consistent Histories 

Formulation 

In both the MDQT and WFR/SP approaches, the classical switching path can be written 
as a time ordered sequence of events, 

i?°(t) = {i?°°,---,i?;^---,i?:j"} (23) 

in which superscript aj denotes the switching outcome at time-step j and Oo = i correspond- 
ing to our choice of the initial quantum state. Two such paths are shown schematically in 
Fig. |1| where we plot the eigenenergy of the occupied state along the path as a function of 
time for an aqueous electron following laser excitation. Changes in the quantum state imply 
that there is a corresponding sudden change in the forces exerted on the classical particles 
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over the course of its evolution. The result is that different sequences of switching events 
will lead to rapidly diverging paths. 

Each path of the quantum subsystem will have associated with it a unique classical path 
characterizing the slight differences in the bath dynamics arising from the slight differences in 
the forces from different quantum paths. These slight differences in the bath dynamics mean 
that bath trajectories will be divergent even on a very short time scale. The result of this 
is that the off-diagonal elements of the quantum density matrix will be rapidly diminished 
as the various bath paths diverge. Thus, coherence in the quantum subsystem is lost due to 
the slight differences in the bath dynamics for each quantum path. We now explore the issue 
of phase coherence between the solvent and solute subsystems. Our approach is based upon 
a novel interpretation of quantum mechanics introduced by Griffiths ||51|, Omnes |]52| , |53 



and Gell-Mann and Hartle [^^ -|56| over roughly the past 10 years. This interpretation of 
quantum mechanics, termed "consistent" or "decoherent" histories, has generated a great 
deal of attention in the field of quantum cosmology [^^-Q. 



Let us consider the evolution of an initial quantum state \i{Ro)), taken as an adiabatic 
eigenstate of the quantum mechanical Hamiltonian for initial bath (or classical) configuration 
Ro, along a switching path, R°'{t). During each time interval, we determine whether or not 
a quantum transition has occurred according to a stochastic switching criteria as discussed 
above, modify the classical dynamics accordingly, and record the outcome of each switching 
attempt. 

We can write the reduced quantum transition probability for a mixed quantum-classical 
system as 

P^^itf) = ( E T,f[R''{t)]T^f[R^t)]p,{Ro)) , 

\{R-{t),Rf^(t)} I „ 

= ( E f/./[i?"(t)]p.(i?o))4[i?^(t)]e'(^[^"l-^[«^l)\ , (24) 

where Pi{Ro) is the probability density of being in the initial state. The average, (■ ■ ■)o, is 
taken over initial configurations, and the sum is over pairs of switching paths 

{R''{t),Rf^{t)\R''iO) = Ro,R^{0) = Ro} (25) 

which start at the initial configuration, Ro, with the quantum state in the initial state and 
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end at any final configuration with the quantum state in state |/) at tf. 

Notice that this definition of the transition probability no longer involves a single path. 
In fact, the transition probability depends upon all other alternative paths which could 
ultimately lead to the final state. When the sequence of switching events for a pair of paths, 
i?"(t) and R^{t) are identical, the contribution from the classical action will vanish and the 
quantum transition probability along a given path is given simply by 

P^f{t) = ( E \U.fmtWp.{Ro)) . (26) 

which is identical to the transition probability we write in the previous section. Applying the 
assumption of Eq. ^ to a switching path implies that the coherences between the quantum 
states are not damped by the bath. This can lead to the overestimation of the non-adiabatic 



transition rates as demonstrated in Ref. 136 



When transitions do occur, we must consider the phase interference between the various 
alternative switching pathways. At short times following a switch, there will be a significant 
contribution from the phase interferences between two alternative paths with similar histories 
and the total transition probability is the sum over the partial transition amplitudes for each 
pathway. However, at longer times following a switch, the action difference between two 
paths will be very large so that the phase interferences will be added destructively, leaving 
a sum over incoherent transition probabilities. In short, there will be considerable phase 
coherence between paths with similar histories but little or no phase coherence between paths 
with dissimilar histories. The timescale over which we must explicitly consider quantum 
interference effects between alternative pathways is the quantum decoherence time. 

Recently we demonstrated that quantum decoherence effects can be consistently incor- 
porated into mixed quantum-classical systems by recognizing that restricting the quantum 
evolution to given classical pathways is equivalent to making a series of quantum measure- 
ments on the total system. |]36| , |48| , ^ The classical path sequence shown in Eq. ^ and illus- 
trated in Fig 1 is an example of a quantum mechanical history. Along this history, the time 
evolution operator for the quantum system can be written equivalently as a time ordered 
sequence of alternating quantum projection operators and unitary evolution operators 

= ?7„P""-i ■ ■ ■ tjiP'^^Uo (27) 
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where the projections at each time interval are members of complete sets representing the 
total set of possible hopping outcomes, 

^P"' = l, (28) 

and 

i,„=,-'/.*:_/-i-<»» (29) 

is the unitary evolution operator for the quantum wavefunction along a segment of the 
switching path. Every time we "roll the dice" to determine the next quantum state along 
the switching path, we need to apply the projection operator P to the quantum wavefunction. 

In general the projection operators in Eq. |28| are any operator acting on the Hilbert space 
of the quantum subsystem. In practice, we have found the Gaussian form to be useful. 

A(Q) = (^)'^'e-"(^-«)'/^ (30) 

where Q is a bounded operator in the Hilbert space of the quantum system with eigenvalues 
{Qj} Later we shall define Q in terms of the coupling between the quantum and classical 
variables. The parameter a serves as the range of values which are projected out by P. 
When acting on the quantum density matrix P acts under the trace operation, 

Pp = (^^y J ciQe-"(^-«)'/2pe-"(Q^Q)V2 

= e-"(«-«^)'/%. (31) 

Thus, when a oo, P projects out all of the adiabatic eigenstate and full coherence between 
the states is retained. On the other hand, taking a = means that we project out only 
the state determined by the random switching event killing off any coherences between the 
eigenstates. 

In our approach, we use the decoherence timescale to set a characteristic time interval 
between subsequent applications of the projection operators. During this interval, transition 
amplitudes are added coherently according to the rules of quantum mechanics. Application 
of the projection operators destroys the coherences and we are left with transition proba- 
bilities which are then added according to the rules of standard probability theory. Within 
our approach, Eq. |4| is rewritten as 
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P.^f{t) = ( E {UC[R-{t)MRo)C^[R'mfn)) 

\B'-{t) I ^ 

where = |/(-R"")) is the final quantum state at the end of path 

When the time intervals between "measurements", 5t„ = t„ — t^-i, are members of a 
Poisson distribution, the probability of maintaining coherence over a given interval, Stn, is 

rSt p-t/r 

7r{5t) = / dt- = 1 - e-^*/", (33) 

JO T 

and T is the characteristic decoherence timescale. Under these conditions, the time evolution 
of the quantum density matrix over the short time step 6t is given as 

Pit + 6t) = (1 - -)(p(t) + i-[HQ,p]) + -Pp. (34) 
r r r 

Taking the limit of 6t ^ yields the master equation for the quantum density matrix 

p = AHq,p\ - -{p-Pp), 

T 

= Cp--{p-Pp). (35) 
r 

Using the Gaussian projector above, 

h = i[HQ,p]- \pij - {Pp\ij). 

= [Cp\ij - -(1 - exp[-a(gi - g,)V4])Pi,. (36) 
r 

The first term contains the Liouvillian of the system, L = i[HQ, ]. Evolution of the quantum 
system under this term alone is unitary and non-dissipative. The second term introduces 
quantum decoherence into the dynamics of the quantum subsystem due to the series of 
projections described above. The coherences originally in the quantum subsystem decay 
due to the series of weak measurements imposed by the environment. 

Both C'[i?"(t)]pi(i?o)C't[i?"(t)] in Eq. |2| and the density matrix in Eq.§B| is a functional 



of only a single classical trajectory as opposed to a sum over pairs of trajectories as in Eq. |24 . 
Where did the other paths go? It turns out that all pairs of paths are still included in these 
equations, except that we periodically "prune" the paths which branch away from the main 



bundle of paths at either the switching times (when we determine a quantum transition to 
have occurred) or when we have determined that a "measurement" has occurred by randomly 
making the density matrix diagonal through application of the projection operators or by 
smoothly diminishing the off diagonal coherence using the master equation. Both methods 
are formally equivalent. 

When the quantum system is coupled linearly in Q to a bath of oscillators, the fluctu- 
ations in the quantum system become proportional to the ability of the bath to dissipate 
energy. Under such conditions for a Markov bath, the relation is [^] 

where 70 is the dissipation constant. One can then relate decoherence in the quantum system 
to fluctuations in the harmonic bath through the fluctuation-dissipation theorem. In essence, 
our decoherence ansatz is entirely consistent with other theories of quantum dissipation. 
In other theories of quantum relaxation, such as the spin-boson model master 



equations such as the Redfield relaxation theory [|65H69|], or Liouville space methods |]T0 



both decoherence and dissipation are treated implicitly through effective equations of motion, 
thus losing the molecular level information about the relaxation dynamics. In our treatment, 
dissipation is included explicitly through the classical molecular dynamics of the condensed 
phase medium, thus retaining a molecular level description. Quantum decoherence is treated 
implicitly in order to avoid summing over alternative pairs of classical paths. 



3. Decoherence Timescales 



Because our formalism relies upon an a priori knowledge of the decoherence time scale 
for the given system, we desire a molecular level theory of the process, allowing access 
to a measure of the quantum decoherence timescale. In developing a theory of quantum 
decoherence we introduced parameters which characterize the decay of coherences due to 
differences in the dynamics of the bath for each possible quantum state. Here we shall 
make use of the projection operator to solidify that concept and determine how one can 
go about obtaining a measure of the coherence time and length scales from realistic mixed 
quantum/classical simulations. 
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Let us consider the effect of coarse graining over a few timesteps. A the initial time to, we 
select the initial state of the total system to be a pure state given by |^'o(-Ro))- Now, evolve 
l^o) to time ti under the time dependent Hamiltonian H{R{t)) where R{t) is a classical 
path and resolve the final state in a basis of eigenstates of H{R{ti)). 

|vI/i(i?0) = EcaJ<(^i))- (38) 
Next, we project this state into coarse grained sets by writing 

Eca,K(it:i)> = E E co^M^iR,)). (39) 

oil Ai aiSAi 

As a result of the projection, there is little or no phase coherence between the various 
eigenstates composing this set. Since each adiabatic eigenstate defines a different adiabatic 
potential which the classical path may follow from this point in time until the next time step, 
each adiabatic state will define a unique coarse grained path or "branch" for the classical 
evolution. Select a branch using some selection criteria and propagate the state under the 
Hamiltonian, H{Ri{t)), where Ri{t) is a classical trajectory along a given branch labeled by 
the sequence of adiabatic potential surfaces followed by the trajectory. We accumulate these 
as an indication of the history of choices made along the path. Again at time t2, resolve 
the new state at t2 in a basis of eigenstates of H{Ri{t2)) and select a new branch for the 
classical dynamics 

E co.MiiRi))-^T. E c.,iCH^2i)). (40) 

aiGAi A2 a2GA2 

Iterating the procedure once again to time ^3 yields the total quantum state of the 
system/bath along a particular trajectory. 

E Ca,|Cn^2l))-E E Ca3lC^^Hi?32l)). (41) 

a2eA2 A3 ctaGAa 

One can begin to see the conceptual underpinnings of the decoherent histories approach 
to quantum mechanics. At each level of coarse graining, the evolution of the classical 
variables branches into different possible histories where each branch indicates evolution 

along a different adiabatic potential surface. The evolution of the total quantum state 
depends very much upon the particular history, {A„, A„_i, • • •}, selected at each time step. 
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In other words, when we allow the classical variables to switch between adiabatic potentials 
at a particular time as the result of a electronic transition, the evolution of the electronic 
system becomes dependent upon when the transitions actually occur. 

Quantum coherence between alternative states is lost as the overlap between the alterna- 
tives decays. We can provide an estimate for the rate at which coherence is lost considering 
the short time evolution of the reduced density matrix along a given path. For example at 
time ^3, the reduced density matrix is given by 



p'^^^^{Rs2,{ts)) = E E E CasCl 

AsA;, aseAs feeA!, 



To estimate the overlap between alternative branches for times t > we assign to the bath 
coordinates at time Gaussian wavepackets centered at the instantaneous phase space 
positions of the bath variable, \G{R, P)). In other words, we write 

I0a?^^'(^32l)) = lC^^^)|G(i?321,P32l)). (42) 

Next, wc make the approximation that we can decouple the evolution of the quantum system 
from the evolution of the bath over a short time period and allow the path to begin to branch 
in to i?32i and i?3'2i- 

p'^'^'{R32i{t>h))= E E 

AsA^ aseAa/JaeA!, 
X lC^^^(i?32l))(0r^n^3'2l)| 

X •'^a3ft(-^321 — -R3'2i; (43) 

J„303(i?321 — i?3'2i; t) is the overlap integral between paths. The decoherence time can thus 
be determined by the timescale at which this vanishes. 

Let us approximate Jas/Ss (t) using Gaussian coherent state wavepackets which follow the 
classical evolution of the alternative paths. Under this approximation, and to lowest order 
in time the overlap decays as (We drop the 321 and 3'21 notation and take R — R321 at 
time ^3 and consider times t > ^3 ). 



Jap{R;t) = exp 



{F"{R) - F'%R)y 



(44) 
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where and F^{R) are the instantaneous quantum forces acting on the classical par- 

ticles as the path bifurcates along the various adiabatic potentials. (Note: that we are 
implicitly summing over all particles in the system in the exponent of this equation.) The 
width of the Gaussian, u, is left arbitrary for the time being. Thus, a rough estimate 
of the quantum decoherence time at a given configuration is given by the decay width of 

j„^(t), mgi 



Finally, averaging over an ensemble of configurations, we obtain the "average coherence 
decay" 



(J(t)) = J dRP{R) exp 



iF''{R)-F''{R)f 



(46) 



4muj 

where P{R) is the probability distribution function for configuration R. 

Such quantities are easily computable in mixed quantum/classical simulations and pro- 
vide a convenient measure for coherence times in simulations. We shall later use this measure 
to assess the choice of decoherence times in simulations of the dynamics of an excess elec- 



tron in ordinary and heavy water. The coherence time for one bath configuration will 
be quite different than the coherence time for another configuration. However, one can envi- 
sion a case in which quantum transitions occur only when the classical variables are in a few 
rare configurations. In this case a global decoherence time would not provide an accurate 
measure of the decoherence timescale. As we shall demonstrate next, even small changes in 
the decoherence timescale can have profound effects on both the quantum and the classical 
evolution. 



4- Coherence Times For Physical Systems 
We can estimate the decay of quantum coherence for the hydrated electron from Eq. 



46| with information available from excited state simulations. For the present example, the 
initial state i is the equilibrium excited state of the hydrated electron, and the final state j 
is the ground state of the electron. For the widths of the frozen Gaussians, we chose 

4muj = GmkT (47) 
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which results from rigorous analysis of the non-adiabatic transition rate between displaced 
harmonic oscillators in the high temperature limit, [^,0 and also allows for direct com- 



parison to the earlier calculations of Neria and Nitzan. p9| , p0| To compute the average de- 
coherence decay we constructed ensemble of excited state configurations by assuming that 
the system was equilibrated at times past 1 ps after being prepared in the excited state, n 
We chose 20 configurations at 25-50 fs intervals from each of 5 individual trajectories for a 
total of 100 configurations. Details of these simulations are given in the Appendix. 

The results of this calculation are shown as in Fig. ^ For the hydrated electron the 
coherence decays in a roughly Gaussian manner, and a Gaussian fit to the decay has a 
variance of ~ 3.1 fs. Another estimate of the decoherence time is found in the area under 
the curve, which for this example is 2.8 fs. This result is in good agreement with previous 



calculations using a different model for the hydrated electron, p9| , p0| and demonstrates a 
posteriori justification for the hypothesis of a ~ 1 fs coherence time in the earlier non- 
equilibrium simulations. We also note that the rapid decoherence of this system provides 
a ^ostenon justification for the short time approximation used to estimate {J{t)). We find 
that the exponential is typically dominated by only 5 to 10 nuclei which are the closest to 
the bulk of the electronic charge density. This makes sense from the stance that the largest 
difference in force between the two surfaces will be for nuclei in positions where the charge 
density, is large on one surface and small on the other. 

We also show the coherence decay curve for an excited electron in D2O computed simi- 
larly. The coherence decay in D2O is qualitatively similar to that in H2O, only for D2O the 
approximate Gaussian decay time is ~ 4.6 fs (versus ;^3.1 fs for H2O) and the area under 
the curve is 4.1 fs (versus 2.8 fs for H2O). 

The longer coherence time in heavy water compared to light water arises predominantly 
from the difference in mass in the choice of the Gaussian width. For classical H2O and D2O, 
the probability of a given nuclear configuration is the same. Static ensemble properties for 
the two fluids should be identical since the ensembles contain identical nuclear configurations 
with equal statistical weights. f\ Since the electronic Hamiltonian for the solvated electron 



^This is true for classical water models. When quantum mechanical effects are properly included, 
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is identical for both heavy and hght water, the static ensemble averaged potential energy 
difference and the difference in Hellmann-Feynman forces on the two surfaces will also be 
identical for the two fluids. Thus, the only differences in the evaluation of the coherence 
decay for the two fluids is through the mass term that enters through the Gaussian width. 
Since the nuclear overlap part of the coherence decay depends on the sum over nuclei, the 
mass change leads to the net slower decay of coherence in D2O versus H2O. In fact, for the 
purposes of evaluating J{t) for the solvated electron in D2O, the H2O simulations would 
suffice. 

The different coherence decay times in the two solvents play a direct role in determin- 
ing the isotope effect on the overall non-adiabatic transition rate. In simplified terms, to 
determine the non-adiabatic transition rate before quantum coherence has decayed, non- 
adiabatic transition amplitudes should be added; after the decoherence interval, memory of 
the complex phases is lost and non-adiabatic transition probabilities should be added. 



5. Effect of Decoherence on Quantum Transition Rates 



One of the primary effects of electronic coherence loss can be observed in the evaluation 
of the switching probabilities. If we approximate the short time evolution of the density 
matrix as 



Pij{t + 5t) ^ exp [1 - exp (-a(Fi - F,)V4)] | 



Pijit), 



(48) 



and substitute this into the non-adiabatic switching probability given above, one obtains 

6t V 



Pij{St) fa exp < 



1 _ exp (-«(/, - /,) 74) 



Jt pii 



(49) 



This predicts that one should see a change over from non-adiabatic dynamics to adiabatic 
dynamics as the coherence timescale becomes very small. Thus, even in this relatively simple 



H2O and D2O have slightly different properties due to differences in the spatial dispersion of the 
H and D nuclei. C. f. G. S. Delbuono, P. J. Rossky, J. Schnitker, J. Chem. Phys. 95, 3728 (1991). 
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model of coherence loss, it is readily apparent that decoherence due to the coupling to an 
external environment generally tends to suppress quantum effects and increases the degree 



of adiabaticity in the quantum subsystem 



As one can see from the last set of equations, one can estimate the effect of coherence 
loss on the quantum transition probability given a coherence timescale and the quantum 
transition amplitudes for a given system. Recall that the origin of this equation is that 
the classical degrees of freedom (the bath) make a series of quantum measurements on the 
system at random time intervals. Thus we can write the non-adiabatic transition probability 
between states i and j as a function of the measurement timescale r as 

2 



n=l 



(50) 



Here, Tj;l^\nAt) is the quantum transition amplitude at the n^^ timestep starting from the 



Tth 



N configuration. In other words, to compute this quantity, we propagate the quantum 
subsystem coherently over a coherence time picked from a distribution of possible coherence 
times. During this interval we sum over amplitudes. At the end of the interval, coherences 
between states are killed off and we sum over probabilities. The outer summation is over 
starting configurations and WN{T,nAt) is the statistical weight for a given configuration to 
have a coherence time of r. These are normalized by Q. The time step. At is taken as the 
shortest timestep used to propagate the equations of motion for the system. For the case of 
the electron in H2O and D2O, this r = 1 fs throughout. 

In Fig. ^, we plot the estimated lifetime of the electron in the excited state as a function of 
coherence timescale for both solvents. The solid lines, labeled "fixed" refer to using a single 
coherence time interval in evaluating Eq.^, i.e. WN{T,nAt)/Q = 9{t — nAt). The curves 
labeled "Poisson" refer to lifetimes computed using a Poisson distribution of coherence time 
intervals. 

wn{t, nAt) IQ = exp(-nAt/r) /r (51) 

Both curves track each other quite well. For a 1 fs coherence time, these lifetimes correspond 
to ~ 550 fs for H2O and ~ 850 fs for D2O. The magnitude of these rates agree reasonably 
with the rates reported in earlier simulations p6|,[7T|,|72| and the roughly 2:1 simulated isotope 
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effect is reproduced as well. While non-adiabatic transitions typically occur from configura- 
tions with higher than average transition probabilities, these configurations are reasonably 
rare that the average transition probability determined here provides a reasonable estimate 
of the non-equilibrium population dynamics. 

We can now use the coherence timescales estimated above to provide an estimate of 
the excited state survival times. For equal coherence times of 1 fs, which was the original 
assumption in the simulations by Schwartz and Rossky, the survival probabilities in Fig.|^ 
predict roughly a 2:1 isotope effect. However, when we use the decoherence times computed 
above, we find that the ~ 50% difference in coherence times yields an estimation of the 
lifetimes which are identical to within 10%. Furthermore, the lifetimes and estimated isotope 
effect agree remarkably with the experimental values reported by Barbara. These results 
are summarized in Table 1. This provides a clear demonstration that quantum decoherence 
plays a direct and key role in the electronic dynamics of this system. 



B. Decoherence in Quantum MD Simulations 

Quantum decoherence can be also be incorporated into non-adiabatic quantum MD sim- 
ulations through either the projection operator method or the master equation theory dis- 
cussed above. In the statistical limit, both methods give identical results. The com- 



putational algorithm given below uses the projection operator approach (Eq.|3^) with the 
coherence time intervals chosen from a Poisson distribution with characteristic time scale 
To = 3.1 fs as derived from our estimates of the coherence time scale for an electron in H2O 
when the excited state is nearly solvated. This time scale was estimated by computing the 
average decay time of the overlap of frozen Gaussian vibrational wavefunctions evolving on 
different adiabatic potential energy surfaces by sampling a large number of nearly equivalent 
excited state configurations. Also note that this method can be used for both the MDQT 
method and the WFR/SP methods since we assume that the treatment of the quantum- 
classical forces is independent of coherence retention. The consistent histories algorithm 
proceeds as follows: 

1 . Determine coherence interval from Poisson distribution of possible intervals with char- 
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acteristic time scale td- 

2. Propagate the quantum wavefunction over this time scale while self consistently evolv- 
ing the classical degrees of freedom. At the end of each dynamical time step h < t^, 
determine the switching path followed by the classical variables using the stationary 



phase algorithm developed by Webster, et al. [p3| , p4| modified such that coherence in 
the quantum wavefunction is maintained over the entire coherence interval. We also 
note that the stationary phase switching path is determined in a piece-wise continuous 
fashion by selecting intermediate quantum states every dynamical time step. This is 
to avoid the computational overhead of computing variationally the stationary phase 
trajectory over relatively long time intervals. So long as td is no longer than a few 
dynamical time steps, this approximation should not be too extreme; however, certain 
pathological cases can be invented in which this approximation does break down. 

3. If either a switch occurs in the time interval or we reach the end of the interval, the 
quantum wavefunction is collapsed using the projection operators discussed above. 

4. Repeat. 

In Fig. ^ we plot the switching times from the excited state to the ground state for 
a total of 23 simulation runs. The starting configurations were generated by performing 
a 35 ps simulation in which the electron was prepared in the ground state and 16 initial 
configurations were chosen whenever the energy difference between the ground state and one 
of the p-like excited states become resonant with the excitation laser (1.7 eV). ||35| , |28|j73 



These starting configurations were identical to configurations used previously by Schwartz 
and Rossky in their work on this system. ||2^,7T,7^ Following the initial excitation, the 



system was allowed to evolve. During this time, the energy gap between the p-states and 
the s-state narrowed to ~ 0.5eV as the excited state was solvated by the surrounding water 
molecules. The simulation continued until a switch from the excited state to the ground 
state was recorded. Immediately after the switch, the energy gap between the occupied 
ground state and the first excited state widened dramatically as the solvent responded to 
the new electronic state. The average and median switching times from these simulations 
are compiled in Table 2. 
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In order to test the sensitivity of the switching times to the choice of sequence of coherence 
intervals, five configurations were "recycled" by using different random number sequences 
for the coherence time intervals to generate different switching paths starting from the same 
initial configuration. For each such pairs of paths, the classical dynamics were identical 
up until the earlier switching time when one of the paths switched from the excited state 
to the ground state. However, since the coherences between the different quantum states 
were killed off at randomly different times for each path, each path sampled a different 
probability distribution function for making a switch at each MD time step. Hence any 
correlations between switching times resulting using the same initial configuration reflect 
correlations in the distribution of coherence time intervals. Given that there is typically a 
100 fs time difference between pairs of data (in one case a 530 fs difference) which is roughly 
1/3 of the average survival time scale of the excited state, we find very little correlation 
between the switching times originating from the same initial configuration. 

As the excited state is solvated by the surrounding water molecules, the energy gap 
between the excited state and the ground state narrows to its equilibrium value. Further- 
more, from first order perturbation theory, we expect that the electronic transition rate to 
be inversely proportional to the magnitude of the energy difference between initial and final 
states. Thus, a simple model for the excited state survival probability can be written as 

subject to the initial condition P{0) = 1. Here, fceq is the non-adiabatic transition rate 
for the solvated excited state and u{t) is the average energy gap normalized to the average 
energy gap of the solvated excited state. Lj{t) is related to the solvation response S{t) by 

i,it) = i-')sm + Mil - s{t)) ^ j^3j 

In Fig. ^, we plot u(t) using the fit to the solvation response function generated by our 
simulations. These fits are nearly identical to those obtained by Schwartz and Rossky. 
?T| , |25|PB| , [7T1 ] At short times, when u!{t) > 1, the non-adiabatic transition rate will be small 



since the energy gap is large. At long times uj{t) ^ 1 as the energy gap relaxes to its 
equilibrium value. Thus, k^q is the exponential decay constant of the excited state population 
once the excited state is fully solvated. 
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The solution of Eq. |5^ is the one parameter family of curves given by 

■.J'^]pi0). (54) 
"'0 E{s) J 

Using a non-linear fitting procedure, we fit our data to this model to obtain an estimate 
of the equilibrium non-adiabatic lifetime of 241 fs with a = 1-707. A plot of our data 
superimposed on this fit is given in Fig. ^. Interestingly enough, our data does not fit this 
simple model as nicely as the data given by Schwartz and Rossky. This suggests that the 
longer coherence times used in this study imparts a non-trivial memory dependency into the 
survival probability which is inadequately captured in Eq. |5^. 

As expected, the lifetimes reported here (Table 2) are consistently shorter than the 
lifetimes reported by Schwartz and Rossky in which a constant 1 fs coherence time scale 
throughout their simulations, thus emphasizing the profound sensitivity of these simulations 
to the coherence time scale. Furthermore, our results are consistent with the estimated 
lifetimes estimated above. 



P(t) = exp 



IV. DISCUSSION 



In this chapter we have briefiy described the results of our work towards a molecular level 
description of quantum relaxation phenomena. Here we have focused exclusively upon the 
role that transient quantum mechanical coherences between the solvent and the solute play in 
the electronic relaxation of an excited solute species. In mixed quantum-classical computer 
simulations, fundamental assumptions about the decay of these coherences produce direct 
manifestations on the computed quantum mechanical transition rates and must be included 
consistently in order to make realistic predictions and comparisons. PB| , |i5[] 

We have focused our attention primarily upon electronically non-adiabatic processes and 
various aspects which must be considered in order to simulate such processes for a condensed 
matter system. In most quantum-chemical calculations, the fundimental assumption is that 
the nuclear dynamics are well within the Born-Oppenheimer approximation and that break- 
downs in this approximation are rare and unusual events. However, in nature the opposite 
is largely true. Non-adiabatic events occur in a number of important chemical processes 
including the process of vision and photosynthesis, predissociation, photodissociation, charge 
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transfer, spin-quenching, and charge transfer reactions of which free radical chemistry is a 
vital component. 

The computational algorithm, based upon the so-called consistent histories interpreta- 
tion of quantum mechanics, provides both the molecular level underpinnings of quantum 
decoherence and the computational means for properly including decoherence effects in 
non-adiabatic quantum-molecular dynamics simulations. According to the rules of ordinary 
quantum mechanics, a quantum system will evolve into a coherent superposition of alter- 
native states. In our decoherence theory, this coherence is dissipated due to the differences 
in the forces exerted on the bath by alternative states involved in the superposition. Thus, 
on short scales, a quantum system in a bath will obey the rules of ordinary quantum me- 
chanics and evolve into a coherent superposition of states, whereas on longer time scales, 
the coherences between states are diminished and the quantum system must be described as 
statistical (i.e., incoherent) mixture of states. As the quantum system interacts continuously 
with the bath, coherences between states are continuously created by non-adiabatic coupling 
and damped by the divergence in the bath dynamics induced by the system-bath coupling. 

This subtle interplay between coupling and decoherence and the subsequent dependency 
of transition rates on the decoherence time scale has profound implications for a variety of 
condensed phase chemical dynamics including: internal conversion and internal vibrational 
energy redistribution (in which the bath is comprised of all the modes of the molecules 
except for the one mode of interest), electronic energy transfer between molecules or different 
parts of the same molecule, and charge transfer reactions including proton and electron 
transfer. In these latter examples, both the condensed phase environment and the internal 
motions of the molecules act as a bath which couples the quantum states together. The 
decay of quantum coherence, which depends upon the frequencies and populations of the 
bath modes coupled to the quantum system will determine the extent to which the non- 
adiabatic coupling can act to allow the chemical reaction to proceed. Changes in the spectral 
density due to isotopic changes in the bath can have a substantial impact on non-adiabatic 
chemical dynamics f^pofj . Furthermore, the decay of quantum coherence can determine the 
adiabaticity for a chemical reaction. 

Perhaps the two major lacunae in our theory of decoherence is the explicit dependency 
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upon an a priori estimate of the coherence time scale and the fact that this time scale 
remains fixed throughout the calculation. As mentioned above, we estimated this time scale 
by computing the average decay time of the overlap of a product of Gaussian coherent states 
evolving on different adiabatic potential energy surfaces. In this estimate, the individual 
widths of the coherent state wavefunctions centered about the initial phase space points of 
the classical nuclei are set to be proportional to the thermal DeBroglie wavelength of each 
nuclei. Although physically realistic for a variety of situations, this does leave the coherence 
length scale (i.e. the widths) as an adjustable parameter. While the effects of changing 
the coherence length scale over a broad range have not been systematically studied, results 
from our previous work demonstrate that the non-adiabatic transition is quite sensitive to 
changes in the coherence time scale and hence will be sensitive to changes in the coherence 
length scales. Furthermore, as the bath explores various regions of the quantum potential 
energy surface, the force differences which give rise to the decay of the quantum coherences 
should vary from one configuration to the next. |HS| Current progress is underway towards 



obtaining both the coherence length and time scale during the course of the non-adiabatic 
simulation ab initio by examining the quantum mechanical fiuctuations of the bath particles 
about their stationary phase paths. 
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Tables 



TABLE I. Summary of experimental and theoretical estimates of the nonradiative lifetime of 
an excess electron in H2O and D2O. The first two rows are experimental results from Gauduel 
and Barbara and the remainder are theoretical estimates. All theoretical numbers reported are 
for relaxation from an equilibrated first excited state of the electron in either solvent. See text for 
details regarding the differences in the theoretical values. 





H2O (fs) 


D20(fs) 


D2O/H2O 


Gauduel et al. (Ref. 0) 


240 


250 


1.04 


Kimura et al. (Ref. ||) 


310 ± 80 


310± 80 


1.0 


Neria and Nitzan (Ref. 


220 


800 


« 4 


Schwartz and Rossky (Ref. PH) 


450 


850 


2 


Ifs Coherence" 


550 


850 


1.55 


Present Estimate* 


310 - 285 


345 - 334 


1.11 - 1.17 



1.0 fs coherence time for both H2O and D2O. 



^ Estimated hfetime ranges correspond to estimated ranges of the coherence time for both 
solvents. H2O: 2.8 -3.1 fs. D2O: 4.1-4.6 fs. (c.f. Fig 0) 



TABLE II. Excited state lifetimes for the aqueous electron. (SR= Schwartz and Rossky, J. 
Chem. Phys. 101 6902 (1994). SBPR = Schwartz, Bittner, Prezhdo, and Rossky, J. Chem. Phys. 
104 4942 (1996)). The coherence times used in each study are listed in parentheses. 





Present (3.1 fs ) 


SR (1 fs) 


SBPR (2.8-3.1 fs)(^) 


Median 


338 fs 


630 




Average 


384 


730 




Equilibrium 


234 


450 


310-270 



(a.) Only Equilibrium lifetimes considered. 
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Figure Captions 



FIG. 1. Energy levels and occupations for two trajectories starting from identical water config- 
urations. The occupied state as a function of time is indicated by + and o for the two trajectories. 

FIG. 2. Overlap function averaged over 100 equilibrated excited state configurations for H2O 
and D2O. 

FIG. 3. Effective Excited State Lifetimes vs. Quantum Decoherence Time. Shown for both H2O 
and D2O are estimates of the excited state lifetimes using the estimates of the quantum coherence 
timcscalcs. As discussed in the text, we estimate the coherence times to be 2.8-3.1 fs for H2O and 
4.1-4.6 fs for D2O. The corresponding lifetimes are 310 - 285 fs and 345-334 fs for H2O and D2O 
respectively. The experimental range of 310ib 80 fs is from Ref. [28]. See text for details. 

FIG. 4. Distribution of switching times following initial excitation for e~/H20. Diamonds (o) 
represent the switching times from the simulations. Superimposed curves are the results of a 
Gaussian fit (B) and a kinetic model (C). 

FIG. 5. Plot of the normalized energy gap between the ground and occupied excited state 
following initial excitation. 
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